Relative contributions

SLA individual for full dataset

RC_SLA<- lmer(SLA100 ~ treatment*elev*mat_treat+S_init_size +(1|block),control = lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e7)), data = gh2)
Anova(RC_SLA)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: SLA100
##                            Chisq Df Pr(>Chisq)    
## treatment                13.1003  1  0.0002952 ***
## elev                     13.6912  1  0.0002155 ***
## mat_treat                 7.3609  2  0.0252121 *  
## S_init_size              38.9269  1    4.4e-10 ***
## treatment:elev            2.3215  1  0.1275961    
## treatment:mat_treat       1.6117  2  0.4467191    
## elev:mat_treat            1.3702  2  0.5040389    
## treatment:elev:mat_treat  0.0353  2  0.9825290    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (RC_SLA)

visreg(RC_SLA, overlay = TRUE, "elev", by="treatment", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

visreg(RC_SLA,overlay = TRUE, "elev", by="mat_treat",type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
      #fill=list(col=grey(c(0.7,0.9))),
      #line=list(col=grey(c(0,0.5))),
       #points=list(cex=1.5,col=grey(c(0.2,0.8)))
      )  

SLA family level

#read in emmeans
LSmeans_SLA <- read.csv("LSmeans_SLA_2.csv", stringsAsFactors=TRUE)

RC_SLA_means<- lm(emmean ~ elevation*treatment*mat_treat, data = LSmeans_SLA)
Anova(RC_SLA_means)
## Anova Table (Type II tests)
## 
## Response: emmean
##                               Sum Sq  Df F value    Pr(>F)    
## elevation                      13411   1  6.7190   0.01095 *  
## treatment                      73602   1 36.8764 2.229e-08 ***
## mat_treat                       5372   2  1.3458   0.26497    
## elevation:treatment             2232   1  1.1185   0.29277    
## elevation:mat_treat             2575   2  0.6450   0.52680    
## treatment:mat_treat              690   2  0.1729   0.84151    
## elevation:treatment:mat_treat     22   2  0.0056   0.99441    
## Residuals                     201587 101                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(RC_SLA_means)

##use this

visreg(RC_SLA_means, overlay = TRUE, "elevation", by="treatment", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

visreg(RC_SLA_means, overlay = FALSE, "elevation", by="mat_treat", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="specific leaf area", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       #fill=list(col=grey(c(0.7,0.9))),
       #line=list(col=grey(c(0,0.5))),
       #points=list(cex=1.5,col=grey(c(0.2,0.8)))
       ) 

Trichomes individual

RC_Tri<- glmer(trichomes ~ mat_treat*elev*treatment+S_init_size +(1|block),family = Gamma(link=log), data = gh2) 
## boundary (singular) fit: see ?isSingular
#had errors but this seems to work when i specify link=log
#singular fit

Anova(RC_Tri) # elev is sig, sig interaction btw mat_treat and elev
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: trichomes
##                             Chisq Df Pr(>Chisq)    
## mat_treat                  9.6468  2   0.008039 ** 
## elev                     111.6048  1  < 2.2e-16 ***
## treatment                  2.2309  1   0.135277    
## S_init_size               15.4506  1  8.469e-05 ***
## mat_treat:elev            10.5228  2   0.005188 ** 
## mat_treat:treatment        1.5863  2   0.452418    
## elev:treatment             0.1045  1   0.746480    
## mat_treat:elev:treatment   0.7294  2   0.694389    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (RC_Tri)

visreg(RC_Tri, "elev", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Trichome density", partial=FALSE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  
## Warning:   Note that you are attempting to plot a 'main effect' in a model that contains an
##   interaction.  This is potentially misleading; you may wish to consider using the 'by'
##   argument.
## Conditions used in construction of plot
## mat_treat: Parent
## treatment: Control
## S_init_size: 0.251658
## block: 9

visreg(RC_Tri, "elev", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  
## Warning:   Note that you are attempting to plot a 'main effect' in a model that contains an
##   interaction.  This is potentially misleading; you may wish to consider using the 'by'
##   argument.
## Conditions used in construction of plot
## mat_treat: Parent
## treatment: Control
## S_init_size: 0.251658
## block: 9

visreg(RC_Tri, "elev",  by="mat_treat", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Trichome density", partial=FALSE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

visreg(RC_Tri, "elev",  by="mat_treat", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="ç", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

Trichomes family level

LSmeans_tri <- read.csv("LSmeans_tri.csv", stringsAsFactors=TRUE)

RC_tri_means<- glm(emmean ~ mat_treat*elev*treatment, data = LSmeans_tri) 
Anova(RC_tri_means)
## Analysis of Deviance Table (Type II tests)
## 
## Response: emmean
##                          LR Chisq Df Pr(>Chisq)    
## mat_treat                  0.5249  2  0.7691807    
## elev                      15.0428  1  0.0001051 ***
## treatment                  0.0234  1  0.8782974    
## mat_treat:elev             2.6707  2  0.2630689    
## mat_treat:treatment        0.3957  2  0.8205087    
## elev:treatment             0.0217  1  0.8827648    
## mat_treat:elev:treatment   0.1180  2  0.9427011    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(RC_tri_means)

##use this

visreg(RC_tri_means, overlay = TRUE, "elev", by="treatment", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

visreg(RC_tri_means, overlay = FALSE, "elev", by="mat_treat", type="conditional", scale = "response", 
       xlab="Elevation (m)", ylab="Trichome density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       #fill=list(col=grey(c(0.7,0.9))),
       #line=list(col=grey(c(0,0.5))),
       #points=list(cex=1.5,col=grey(c(0.2,0.8)))
       ) 

Plasticity overview

SLA

#read in emmeans. using lsmeans because its easier to see the graph
LSmeans_SLA <- read.csv("LSmeans_SLA_2.csv", stringsAsFactors=TRUE)

SLA <- lm(emmean ~ treatment*mat_treat, data = LSmeans_SLA)
Anova (SLA)
## Anova Table (Type II tests)
## 
## Response: emmean
##                     Sum Sq  Df F value    Pr(>F)    
## treatment            74246   1 36.1403 2.587e-08 ***
## mat_treat             5233   2  1.2736    0.2840    
## treatment:mat_treat    652   2  0.1586    0.8535    
## Residuals           219818 107                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##use this

visreg(SLA, overlay = TRUE, "mat_treat", by="treatment", type="conditional", scale = "response", 
       xlab="Maternal environment", ylab="Specific leaf are", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

#plasticity figure
ggplot(
  #subset(
    LSmeans_SLA,
   # treatment == "Herbivorized"),
  aes(x=treatment, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Specific Leaf Area") +
  xlab("Treatment")
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Removed 1 rows containing non-finite values (stat_summary).
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 38 rows containing missing values (geom_segment).

#plasticity figure
ggplot(
  subset(
    LSmeans_SLA,
   treatment == "Herbivorized"),
  aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Specific Leaf Area") +
  xlab("Maternal enviornment")
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 57 rows containing missing values (geom_segment).

#plasticity figure
ggplot(
  subset(
    LSmeans_SLA,
   treatment == "Control"),
  aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Specific Leaf Area") +
  xlab("Maternal enviornment")
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 56 rows containing missing values (geom_segment).

Trichome Density

#read in emmeans. using lsmeans because its easier to see the graph
LSmeans_tri <- read.csv("LSmeans_tri.csv", stringsAsFactors=TRUE)

TRI <- glm(emmean ~ treatment*mat_treat, data = LSmeans_tri)
Anova (TRI)
## Analysis of Deviance Table (Type II tests)
## 
## Response: emmean
##                     LR Chisq Df Pr(>Chisq)
## treatment            0.00821  1     0.9278
## mat_treat            0.43504  2     0.8045
## treatment:mat_treat  0.28404  2     0.8676
##use this

visreg(TRI, overlay = TRUE, "mat_treat", by="treatment", type="conditional", scale = "response", 
       xlab="Maternal environment", ylab="Trichome Density", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col=grey(c(0.7,0.9))),
       line=list(col=grey(c(0,0.5))),
       points=list(cex=1.5,col=grey(c(0.2,0.8))))  

#plasticity figure
ggplot(
  #subset(
    LSmeans_tri,
   # treatment == "Herbivorized"),
  aes(x=treatment, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Trichome Density") +
  xlab("Treatment")
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Removed 1 rows containing non-finite values (stat_summary).
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 38 rows containing missing values (geom_segment).

#plasticity figure
ggplot(
  subset(
    LSmeans_tri,
   treatment == "Herbivorized"),
  aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Trichome Density") +
  xlab("Maternal enviornment")
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 57 rows containing missing values (geom_segment).

#plasticity figure
ggplot(
  subset(
    LSmeans_tri,
   treatment == "Control"),
  aes(x=mat_treat, y=emmean, color = genotype)) + # Change fill to color
  theme_classic(base_size = 22) + 
  #geom_point() + 
  stat_summary(fun=mean, position = "dodge") + 
  #stat_summary(
  #  geom="errorbar", 
  #  fun.data= mean_cl_boot,
  #  width = 0.1, size = 0.2, col = "grey57"
  #  ) + 
  # Lines by species using grouping
  #scale_color_viridis(discrete = FALSE) +
  #scale_fill_manual(values = cbp2)+
  stat_summary(aes(group = genotype), geom = "line", size=1.5, fun = mean) +
  ylab("Trichome Density") +
  xlab("Maternal enviornment")
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Warning: Width not defined. Set with `position_dodge(width = ?)`
## Warning: Removed 56 rows containing missing values (geom_segment).

Fitness models: Fitness ~ garden x treatment x elev {.tabset)

Probability of reproduction

Individual

mod_flower<- glmer(flowered~ elev*treatment*mat_treat+ S_init_size + (1|block)+(1|genotype), data = gh2, control=glmerControl(optimizer="optimx", tolPwrss=1e-3, optCtrl=list(method="nlminb")), family=binomial)
## Loading required namespace: optimx
Anova(mod_flower) # all seperate factors are sig, sig interaction btw mat_treat and elev
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: flowered
##                            Chisq Df Pr(>Chisq)    
## elev                      1.6336  1   0.201202    
## treatment                10.2749  1   0.001348 ** 
## mat_treat                10.1867  2   0.006138 ** 
## S_init_size              71.0068  1  < 2.2e-16 ***
## elev:treatment            1.3670  1   0.242324    
## elev:mat_treat            5.4014  2   0.067158 .  
## treatment:mat_treat       0.1544  2   0.925696    
## elev:treatment:mat_treat  4.8511  2   0.088428 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot (mod_flower) #plot is strange, not sure if ok fit

flower <-predictorEffect("elev",  partial.residuals=TRUE, mod_flower)
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vector
plot(flower, lwd=2,xlab="Source Elevation (KM)", ylab="Probability of reproduction", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), 
          partial.residuals=list(smooth=TRUE, pch=19, col="black"), ylim=c(0,1))
## Warning in plot.eff(flower, lwd = 2, xlab = "Source Elevation (KM)", ylab
## = "Probability of reproduction", : partial residuals are not displayed in a
## multiline plot

#no parents
mod_pr1<- glmer(flowered~ elev*treatment*mat_treat+ S_init_size+(1|block)+(1|genotype), data =subset(gh2, mat_treat!="Parent"), control=glmerControl(optimizer="optimx", tolPwrss=1e-3,optCtrl=list(method="nlminb")), family=binomial)

Anova(mod_pr1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: flowered
##                            Chisq Df Pr(>Chisq)    
## elev                      2.7394  1   0.097903 .  
## treatment                10.8137  1   0.001008 ** 
## mat_treat                 0.0000  1   0.999811    
## S_init_size              47.3706  1  5.876e-12 ***
## elev:treatment            0.0731  1   0.786891    
## elev:mat_treat            5.1058  1   0.023846 *  
## treatment:mat_treat       0.0705  1   0.790671    
## elev:treatment:mat_treat  2.2386  1   0.134606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
flowerP <-predictorEffect("elev",  partial.residuals=TRUE, mod_pr1)
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vector
plot(flowerP, lwd=2,xlab="Source Elevation (KM)", ylab="Probability of reproduction", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,1))
## Warning in plot.eff(flowerP, lwd = 2, xlab = "Source Elevation (KM)", ylab
## = "Probability of reproduction", : partial residuals are not displayed in a
## multiline plot

family level (not working)

Probability of survival

individual

mod_surv<- glmer(survival~elev*treatment*mat_treat + (1|block)+(1|genotype), data = gh2, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e7)), family=binomial)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
Anova(mod_surv) #only treatment is sig
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: survival
##                           Chisq Df Pr(>Chisq)   
## elev                     0.3383  1    0.56083   
## treatment                9.2003  1    0.00242 **
## mat_treat                1.8540  2    0.39574   
## elev:treatment           0.4366  1    0.50878   
## elev:mat_treat           0.5099  2    0.77497   
## treatment:mat_treat      0.0307  2    0.98477   
## elev:treatment:mat_treat 0.0002  2    0.99992   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mod_surv)

survived <-predictorEffect("elev",  partial.residuals=FALSE, mod_surv)
plot(survived, lwd=2,xlab="Elevation of origin", ylab="Fitness (survival)", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), 
          partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,1))

family level not working

Fecundity

individual

library(glmmTMB)
## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.7.21
## Current TMB version is 1.7.22
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
## 
## Attaching package: 'glmmTMB'
## The following object is masked from 'package:gamlss':
## 
##     refit
hurdle <- glmmTMB (Mature_silique_number ~ treatment*mat_treat*elev+S_init_size  + (1| block)+(1|genotype), zi=~ treatment*mat_treat*elev+S_init_size  + (1| block)+(1|genotype),data = gh2 ,family=truncated_nbinom2)

diagnose(hurdle)
## Unusually large coefficients (|x|>10):
## 
##                           zi~(Intercept) 
##                                -12.75103 
##                         zi~mat_treatherb 
##                                 11.52209 
##                       zi~mat_treatParent 
##                                 10.04355 
##   zi~treatmentHerbivorized:mat_treatherb 
##                                -11.35602 
## zi~treatmentHerbivorized:mat_treatParent 
##                                -11.80556 
## 
## Large negative coefficients in zi (log-odds of zero-inflation),
## dispersion, or random effects (log-standard deviations) suggest
## unnecessary components (converging to zero on the constrained scale);
## large negative and/or positive components in binomial or Poisson
## conditional parameters suggest (quasi-)complete separation
## 
## 
## Unusually large Z-statistics (|x|>5):
## 
##                      treatmentHerbivorized 
##                                  -5.370225 
##      treatmentHerbivorized:mat_treatParent 
##                                  -7.215594 
##                 treatmentHerbivorized:elev 
##                                  24.781566 
## treatmentHerbivorized:mat_treatParent:elev 
##                                  11.189258 
##              zi~treatmentHerbivorized:elev 
##                                  -5.921320 
##                            theta_1|block.1 
##                                  -6.004913 
## 
## Large Z-statistics (estimate/std err) suggest a failure of the Wald
## approximation - often also associated with parameters that are at or
## near the edge of their range (e.g. random-effects standard deviations
## approaching 0).  While the Wald p-values and standard errors listed in
## summary() are unreliable, profile confidence intervals (see
## ?confint.glmmTMB) and likelihood ratio test p-values derived by
## comparing models (e.g. ?drop1) may still be OK.  (Note that the LRT is
## conservative when the null value is on the boundary, e.g. a variance or
## zero-inflation value of 0 (Self and Liang 1987; Stram and Lee 1994;
## Goldman and Whelan 2000); in simple cases the p-value is approximately
## twice as large as it should be.)
Anova(hurdle,type="III", component="zi")
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: Mature_silique_number
##                            Chisq Df Pr(>Chisq)    
## (Intercept)               8.4004  1   0.003751 ** 
## treatment                 0.2982  1   0.585032    
## mat_treat                 8.3993  2   0.015001 *  
## elev                      9.2922  1   0.002301 ** 
## S_init_size              56.2618  1  6.344e-14 ***
## treatment:mat_treat       2.9516  2   0.228596    
## treatment:elev            0.0285  1   0.865890    
## mat_treat:elev            8.4809  2   0.014401 *  
## treatment:mat_treat:elev  3.0445  2   0.218217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(hurdle) #zero inflated model is significant, but not sure if i am plotting it?
##  Family: truncated_nbinom2  ( log )
## Formula:          
## Mature_silique_number ~ treatment * mat_treat * elev + S_init_size +  
##     (1 | block) + (1 | genotype)
## Zero inflation:                         
## ~treatment * mat_treat * elev + S_init_size + (1 | block) + (1 |     genotype)
## Data: gh2
## 
##      AIC      BIC   logLik deviance df.resid 
##   1961.9   2108.2   -949.9   1899.9      799 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  block    (Intercept) 0.03159  0.1777  
##  genotype (Intercept) 0.21906  0.4680  
## Number of obs: 830, groups:  block, 12; genotype, 19
## 
## Zero-inflation model:
##  Groups   Name        Variance Std.Dev.
##  block    (Intercept) 0.9225   0.9605  
##  genotype (Intercept) 1.1907   1.0912  
## Number of obs: 830, groups:  block, 12; genotype, 19
## 
## Dispersion parameter for truncated_nbinom2 family (): 4.47 
## 
## Conditional model:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                 5.50261    1.77175   3.106   0.0019
## treatmentHerbivorized                      -0.47436    2.54744  -0.186   0.8523
## mat_treatherb                              -1.83915    1.48197  -1.241   0.2146
## mat_treatParent                            -0.84673    1.60535  -0.527   0.5979
## elev                                       -1.09470    0.58242  -1.880   0.0602
## S_init_size                                 0.09895    0.07036   1.406   0.1597
## treatmentHerbivorized:mat_treatherb         1.45182    3.64446   0.398   0.6904
## treatmentHerbivorized:mat_treatParent      -0.50276    3.62770  -0.139   0.8898
## treatmentHerbivorized:elev                  0.03551    0.88005   0.040   0.9678
## mat_treatherb:elev                          0.64279    0.49719   1.293   0.1961
## mat_treatParent:elev                        0.31005    0.53393   0.581   0.5615
## treatmentHerbivorized:mat_treatherb:elev   -0.47540    1.26443  -0.376   0.7069
## treatmentHerbivorized:mat_treatParent:elev  0.11147    1.24729   0.089   0.9288
##                                              
## (Intercept)                                **
## treatmentHerbivorized                        
## mat_treatherb                                
## mat_treatParent                              
## elev                                       . 
## S_init_size                                  
## treatmentHerbivorized:mat_treatherb          
## treatmentHerbivorized:mat_treatParent        
## treatmentHerbivorized:elev                   
## mat_treatherb:elev                           
## mat_treatParent:elev                         
## treatmentHerbivorized:mat_treatherb:elev     
## treatmentHerbivorized:mat_treatParent:elev   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                                            Estimate Std. Error z value Pr(>|z|)
## (Intercept)                                -12.7510     4.3994  -2.898  0.00375
## treatmentHerbivorized                        3.0553     5.5953   0.546  0.58503
## mat_treatherb                               11.5221     4.1898   2.750  0.00596
## mat_treatParent                             10.0436     4.2318   2.373  0.01763
## elev                                         4.3293     1.4202   3.048  0.00230
## S_init_size                                 -1.2138     0.1618  -7.501 6.34e-14
## treatmentHerbivorized:mat_treatherb        -11.3560     7.6944  -1.476  0.13997
## treatmentHerbivorized:mat_treatParent      -11.8056     7.8773  -1.499  0.13395
## treatmentHerbivorized:elev                  -0.3140     1.8594  -0.169  0.86589
## mat_treatherb:elev                          -3.8702     1.3659  -2.833  0.00461
## mat_treatParent:elev                        -3.0781     1.3756  -2.238  0.02524
## treatmentHerbivorized:mat_treatherb:elev     3.9197     2.5866   1.515  0.12967
## treatmentHerbivorized:mat_treatParent:elev   3.9602     2.6349   1.503  0.13284
##                                               
## (Intercept)                                ** 
## treatmentHerbivorized                         
## mat_treatherb                              ** 
## mat_treatParent                            *  
## elev                                       ** 
## S_init_size                                ***
## treatmentHerbivorized:mat_treatherb           
## treatmentHerbivorized:mat_treatParent         
## treatmentHerbivorized:elev                    
## mat_treatherb:elev                         ** 
## mat_treatParent:elev                       *  
## treatmentHerbivorized:mat_treatherb:elev      
## treatmentHerbivorized:mat_treatParent:elev    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fruit <-predictorEffect("elev",  partial.residuals=TRUE, hurdle,component="zi")
## Warning in Effect.glmmTMB(ans, mod, x.var = 1, xlevels = xlevels, ...):
## overriding variance function for effects: computed variances may be incorrect
## Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
## predictor S_init_size is a one-column matrix that was converted to a vector
plot(fruit, lwd=2,xlab="Source Elevation (KM)", ylab="Fruit", pch=19, type="response",lines=list(multiline=TRUE, lty=2:1, col="black"), partial.residuals=list(smooth=FALSE, pch=19, col="black"), ylim=c(0,25))
## Warning in plot.eff(fruit, lwd = 2, xlab = "Source Elevation (KM)", ylab =
## "Fruit", : partial residuals are not displayed in a multiline plot

Family level

LAR

individual

#have to curate the dataset to focus on herb treatment and remove NAs
#LAR 3 has the most damage, LAR 13 is the most recent

LAR_data <- subset(gh2, treatment == "Herbivorized")

LAR_3 <- drop_na(LAR_data,LAR_prop3) #this removes the rows without LAR values
LAR_3 <- dplyr::select(LAR_3, LAR_prop3, elev, genotype, treatment, mat_treat, exp_ID, S_init_size, block)

gamlss_mod1<- gamlss (formula=LAR_prop3~elev*mat_treat+ S_init_size+ random(block)+random(genotype), family=BE, data= LAR_3)
## GAMLSS-RS iteration 1: Global Deviance = -317.6556 
## GAMLSS-RS iteration 2: Global Deviance = -333.6844 
## GAMLSS-RS iteration 3: Global Deviance = -334.1667 
## GAMLSS-RS iteration 4: Global Deviance = -334.1874 
## GAMLSS-RS iteration 5: Global Deviance = -334.1902 
## GAMLSS-RS iteration 6: Global Deviance = -334.191
summary(gamlss_mod1)
## ******************************************************************
## Family:  c("BE", "Beta") 
## 
## Call:  gamlss(formula = LAR_prop3 ~ elev * mat_treat + S_init_size +  
##     random(block) + random(genotype), family = BE,      data = LAR_3) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)   
## (Intercept)          -1.270656   0.774938  -1.640  0.10192   
## elev                  0.326973   0.249718   1.309  0.19122   
## mat_treatherb        -0.038300   1.110804  -0.034  0.97251   
## mat_treatParent       0.066532   1.096706   0.061  0.95166   
## S_init_size          -0.124692   0.039976  -3.119  0.00196 **
## elev:mat_treatherb   -0.024684   0.358378  -0.069  0.94512   
## elev:mat_treatParent -0.004827   0.353084  -0.014  0.98910   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  logit
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.58078    0.04473  -12.98   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  399 
## Degrees of Freedom for the fit:  29.87062
##       Residual Deg. of Freedom:  369.1294 
##                       at cycle:  6 
##  
## Global Deviance:     -334.191 
##             AIC:     -274.4498 
##             SBC:     -155.297 
## ******************************************************************
plot(gamlss_mod1)

## ******************************************************************
##        Summary of the Quantile Residuals
##                            mean   =  -0.006089556 
##                        variance   =  0.9862509 
##                coef. of skewness  =  0.2791661 
##                coef. of kurtosis  =  5.366551 
## Filliben correlation coefficient  =  0.9871931 
## ******************************************************************
LAR_13 <- drop_na(LAR_data,LAR_prop13) #this removes the rows without LAR values
LAR_13 <- dplyr::select(LAR_13, LAR_prop13, elev, genotype, treatment, mat_treat, exp_ID, S_init_size, block)

gamlss_mod2<- gamlss (formula=LAR_prop13~elev*mat_treat+ random(block)+random(genotype), family=BEZI, data= LAR_13)
## GAMLSS-RS iteration 1: Global Deviance = -608.1084 
## GAMLSS-RS iteration 2: Global Deviance = -745.3908 
## GAMLSS-RS iteration 3: Global Deviance = -864.3819 
## GAMLSS-RS iteration 4: Global Deviance = -963.2275 
## GAMLSS-RS iteration 5: Global Deviance = -1025.301 
## GAMLSS-RS iteration 6: Global Deviance = -1049.247 
## GAMLSS-RS iteration 7: Global Deviance = -1054.845 
## GAMLSS-RS iteration 8: Global Deviance = -1055.889 
## GAMLSS-RS iteration 9: Global Deviance = -1056.105 
## GAMLSS-RS iteration 10: Global Deviance = -1056.16 
## GAMLSS-RS iteration 11: Global Deviance = -1056.175 
## GAMLSS-RS iteration 12: Global Deviance = -1056.18 
## GAMLSS-RS iteration 13: Global Deviance = -1056.181 
## GAMLSS-RS iteration 14: Global Deviance = -1056.182
summary(gamlss_mod2)
## ******************************************************************
## Family:  c("BEZI", "Zero Inflated Beta") 
## 
## Call:  gamlss(formula = LAR_prop13 ~ elev * mat_treat + random(block) +  
##     random(genotype), family = BEZI, data = LAR_13) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)   
## (Intercept)           -2.6303     0.9162  -2.871   0.0044 **
## elev                  -0.1151     0.2955  -0.390   0.6972   
## mat_treatherb         -0.9045     1.2696  -0.712   0.4768   
## mat_treatParent        0.7480     1.2484   0.599   0.5496   
## elev:mat_treatherb     0.2692     0.4091   0.658   0.5111   
## elev:mat_treatParent  -0.2596     0.4023  -0.645   0.5193   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   3.3235     0.0894   37.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Nu link function:  logit 
## Nu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.6009     0.2262   -11.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## NOTE: Additive smoothing terms exist in the formulas: 
##  i) Std. Error for smoothers are for the linear effect only. 
## ii) Std. Error for the linear terms maybe are not accurate. 
## ------------------------------------------------------------------
## No. of observations in the fit:  304 
## Degrees of Freedom for the fit:  22.51564
##       Residual Deg. of Freedom:  281.4844 
##                       at cycle:  14 
##  
## Global Deviance:     -1056.182 
##             AIC:     -1011.151 
##             SBC:     -927.4594 
## ******************************************************************
plot(gamlss_mod2)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.196609 
##                        variance   =  0.6954236 
##                coef. of skewness  =  -0.1419971 
##                coef. of kurtosis  =  2.991946 
## Filliben correlation coefficient  =  0.9931949 
## ******************************************************************
visreg(gamlss_mod1, 'elev', by= "mat_treat", overlay = FALSE, type="conditional", 
       #scale = "response", 
       xlab="Elevation (m)", ylab="Percentage leaf area herbivorized", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col="blue"),
       line=list(col=grey(c(0,0.8))),
       points=list(cex=1.5,col=grey(c(0,0.8)))) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)

visreg(gamlss_mod2, 'elev', by= "mat_treat", overlay = FALSE, type="conditional", 
       #scale = "response", 
       xlab="Elevation (m)", ylab="Percentage leaf area herbivorized", partial=TRUE, cex.lab = 1.5,cex.axis = 1.5,
       fill=list(col="blue"),
       line=list(col=grey(c(0,0.8))),
       points=list(cex=1.5,col=grey(c(0,0.8)))) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## new prediction 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6) 
## New way of prediction in random()  (starting from GAMLSS version 5.0-6)